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Abstract. Shaped laser pulses are a powerful tool to induce population transfer 
between electronic molecular states, and time-dependent perturbation theory is 
suitable for a description of such a transfer in weak external fields. The application of 
perturbation theory in numerical simulations of field matter interactions can lead to 
divergences. In a recent paper [K. Renziehausen et. al. J. Phys. B: At. Mol. Opt. 
Phys., 42:195402, 2009] we explained that the arising error in the norm of the wave 
function can be split into two parts. The first part is related to numerical errors caused 
by the discretisation of time that is required in the simulation and can be suppressed for 
a sufficiently small time step or abolished for an adequate numerical implementation 
of perturbation theory. The second part may cause divergences and is associated 
with the perturbative expansion order. We presented numerical evidence without any 
analytical proof. Here we are focussing on the derivation of analytical expressions to 
interpret the behavior of what we have called in the above mentioned paper 'simple 
algorithm'. The derivation of analytical expressions for the interpretation of what we 
have called in the above mentioned paper 'improved algorithm' are given in another 
paper [K. Renziehausen. arXiv, in preparation, 2012]. Moreover, we introduce here 
a gedankenexperiment to illustrate the influence of the different orders on the field- 
molecule interaction. 
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1. Introduction 

In a recently published work we discussed in how far numerical algorithms for the 
implementation of perturbation theory are applicable to ultra-short laser pulse molecule 
interactions. This is an important issue because time-dependent perturbation theory 
is most commonly employed for the investigation of the interaction of atoms [2] and 
molecules [3] with electromagnetic fields. However, the disadvantage of the application 
of time-dependent perturbation theory is that it is generally not norm-conserving. 
Although non-perturbative and norm-conserving algorithms can be used to solve the 
time-dependent Schrodinger equation, the advantage of perturbative methods is: They 
allow for a clear decomposition of multi-photon processes into contributions which stem 
from different orders. As an interesting example we mention time-resolved four-wave 
mixing spectroscopy, where one determines the third-order polarization associated with 
the signal emitted into a given direction It is possible to analyse four- wave mixing 
experiments with non-perturbative methods [5], El [7] . Nevertheless, within perturbation 
theory it is possible to differentiate contributions form higher order interactions [H |9l [10] . 

In the analysis presented in [1] we used a model of two molecular electronic states 
(ground state (|0)) and exited state (|1))). In this model initially only one of these 
two states is populated. Such conditions are realized, e.g., in pump /shaped- dump 
experiments which have recently been reported [HI [121 113] • In such processes where 
high intensity laser pulses interact with molecules, it is important to understand how the 
results of perturbation theory converge to the exact results. Although this convergence 
behaviour, depending in detail on the chosen numerical parameters, on the molecular 
system and on the electric fields, cannot easily be discribed quantitatively, it is still 
possible to analyse qualitatively the general trends how the deviations of the norm from 
unity depend on the chosen parameters and the physical situation. In our former analysis 
we used a simple and an improved algorithm, and we stated that the norm deviations 
can be decomposed into two parts, which we named the stationary orders and the 
oscillatory orders. These different parts differ in their behaviour when the following 
parameters are changed: The time step for the discretisation of time in the numerical 
algorithm, the order of perturbation in which the wave functions are calculated, the 
shape of the potentials of the two electronic states |0) and |1) and the electric field. 
In short the stationary orders occur only for the simple algorithm; they depend on 
the electric field and the time step; they converge to zero for the limit that the time 
step goes to zero, but they do not depend on the shape of the potentials. Thus, we 
concluded that norm deviations caused by the stationary orders are purely numerical 
and they are not related to norm deviations that are caused by the discarding of higher 
order interactions. In contrast to the norm deviations caused by the stationary orders 
the norm deviations caused by the oscillatory orders occur for both algorithms, they do 
not depend in leading order on the time step and thus they do not vanish in the limit 
of a vanishing time step. Moreover, they depend strongly on the chosen perturbation 
order for the wave function, the shape of the potentials of the two electronic states |0) 
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and |1), and on the electric field. This clearly indicates that the errors are related to 
the perturbative expansion of the wave function. The numerical results presented in 
[T] were interpreted with the help of analytical expressions which, however, were given 
without proofs (in particular the equations (13) and (14) of [1] are not proven there). 
It is the purpose of the present paper to fill this gap for those analytical expressions, 
which are related to the simple algorithm. The calculations for the improved algorithm 
will be presented in another paper [T^ . Besides the mathematical analysis of the norm 
deviations of the perturbative wave functions new numerical results are presented, and 
an interpretation in terms of a scattering gedankenexperiment are given. 

The paper is organized as follows: We describe in Sec. [2] the structure of the 
discussed Hamiltonian, summarize the basis of perturbation theory and show how the 
simple algorithm introduced in [l] using perturbation theory for the calculation of the 
wave function can be derived. Section [3] contains an analytical analysis of the wave 
function calculated with the simple algorithm, and in Sec. H] an analytical analysis of 
the norm of this wave function is given (in this section we derive the above mentioned 
equations (13) and (14) of jl] for the simple algorithm). Then, in Sec. [5] we summarize 
the interpretation of the analytical results of Sec. H] and provide the above mentioned 
gedankenexperiment. The paper is finished by a summary in Sec. [61 



2. Theory 

2.1. Hamiltonian 

As mentioned in the introduction, we investigate the interaction of an ultrashort laser 
pulse with a molecule in a model where we consider two electronic states |1) and |0). 
The nuclear degrees of freedom are represented by a single coordinate R. The total 
Hamiltonian H{R,t) consists of the system Hamiltonian Ho{R), and the field-matter 
interaction term W{t) 

H{R,t)=Ho{R) + W{t) 

f + Vi{R) \ f -fiE{t) 

f + VoiR) / \ -fiE{t) 



(1) 



where T is the kinetic energy operator, and where Vj{R), j G {0, 1} are the potentials in 
the electronic states. The perturbation consists of the dipole interaction with the electric 
field E(t) of the laser pulse and the projection /x of the transition dipole moment on the 
laser polarization vector. We take into account the Condon approximation and neglect 
the dependence of the transition dipole-moment on the nuclear coordinates. Moreover, 
dipole-coupling within a single electronic state is not regarded. As we analyse a system 
with two electronic states, we have to work with a two-component nuclear wave function 
"^{Rjt) and the time- dependent Schrodinger equation reads (in atomic units): 
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Due to the structure of the perturbation operator W{t) for even powers of W{t) with 
?7 G N the following equation holds: 

W{tf'' = fi^'^Eitf^ ^ J ° j = ^^^E{tfn (3) 

In order to emphasize this point, we will denote even powers of the perturbation operator 
W{t) without an operator head as W{tY'^ in the following. According to ^ for odd 
powers of the perturbation operator W{t) we have: 

iy(t)2^-i = (4) 

=^ [W{tf^^\H,{R)] = W{tf^^-'\W{t),H,{R)] (5) 

Thus odd powers written in the form W^^'^ are proportional to the perturbation 
operator W{t). As initial condition we fix [T| [T5]: 

«(«.«^«^(*'<^„'^*°') (6) 

In the following calculations, for clarity, we suppress in the notation dependencies on 
the vibrational coordinate R for all quantities. 

2.2. Perturbation theory 

The starting point of time-dependent perturbation theory is the integral equation for 
the wave function [16], 



t 







where we assumed that the interaction starts at time to = 0. In perturbation theory 
the wave function "^{t) is expanded in orders of the interaction-operator, and an 
approximative wave function "^{t, k) is obtained which contains all terms up to order k. 

The wave function in first order \E^(t, 1) is obtained by substituting the exact wave 
function ^(t') appearing under the integral by the initial function (0^^-order wave 
function) evolving in time with the system propagator (\E'(t, 0) = e^*^''*\E'(0)): 

i 

^(t, 1) = e-'^o* ^(0) -ijdt' e-'^o(*-*') W{t') ^(t', 0) (8) 



By iterating ([8]) we obtain higher-order corrections as: 

t 

^(t, k) = e-'^o* ^(0) -t J dt' e-'^°(*-*') W{t') k-1) (9) 







We can use ([9]) as a basis to devise a numerical algorithm for the calculation of 
perturbative wave functions [17]. For this aim we discretise the time t into time steps 
At yielding a time-grid where the times t„ are defined as t„ = n At with whole-number 
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values for n. As a next step we set up an iteration scheme, where k) is calculated 

from '^{tn, k). 

Therefore the integral in ([9]) is divided into a first integral in the limits from t' = to 
t' = in, and a second one reaching from t' = t„ and t' = t„ + At, yielding: 

tn+At 

-I J rft'e-*^°(*"+^*-*') ly(t') (10) 

in 

This expression ( ITO|) for the wave function can be interpreted easily: The first term 
represents an unperturbed time-evolution of the system (during the time-interval 
[tn, tn + At]), whereas the second term stands for the possibility that during this interval 
at least one interaction takes placejjl With these considerations a numerical scheme for 
the evaluation of the wave functions can be developed, which we call the simple algorithm 
(S). This algorithm is constructed by a replacement of the integral in (fTOl) by a single 
term at the time t' = t„+i, what leads to wave functions \E'5'(t„, At; fc) depending on 
tn, At and k. Introducing the abbreviating notations ^'(n, /c) := \E's'(t„, At; /c) and 
W{n) := W{tn), we get the following equation for the simple algorithm [1]: 

^s{n + l,k) = e-^^o^* ^s{n, k) - lAt W{n + 1) ^s{n + 1, A; - 1) (11) 

with the start conditions \l/5(0, k) = \l/(0). 

As a last purpose in this chapter we explicate how the short-time propagator 
^-zHoAt^ which appears in f|TT]) and moves the wave function ^^(r;,,^) over a time step 
At, can be executed numerically: This is customarily done by the split operator method 
of Feit and Fleck [18], where a grid for the spatial coordinate R is used. Being a one 
step method correct in second order in the time step At, the application of the split- 
operator method in the simple algorithm 5* does not diminish the order in the time-step 
At, in which the simple algorithm is correct. The reason for this is as we will see 
in the following Sec. [3] that the simple algorithm is a one step method which applies 
perturbation theory correctly only in first order in At. 

— * 

3. Error analysis of the wave functions \l/s'(n, k) 

In this analysis of the wave functions "^si^^k), first we state a closed form for them, 
and then we show that the simple algorithm is a one step method correct in first order 
in At. 

Before we start our analytic analysis of the wave functions "^si^jk), we have to 
introduce notations that are important for the following. 

First, we define the sequence of non-commuting operators Aj when we use the product 

I Here wc note that of Rcf. [1], Eqn. (9), which corresponds to (fTO]) m this paper, is misleading, 
because there interactions taking place in the small time interval [i„,t'] were discarted. 
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symbol Y\: 



Ha 



AqAi ■■■An 



(12) 



In the subsequent calculations combinatorial arguments are important. In particular, we 
will analyse combinatorial problems, where we have to calculate sums over all possible 
combinations with repetition. In these combinations m elements are taken out of a set 
that contains n elements, and the sequence in which the elements are chosen has no 
relevance. We specify a particular combination with repetition by a vector i7("'"*) that 
has n components, which are natural numbers or zero. The j-th component z/j"'"^'* of 
such a vector equals the number of cases how often the j-th element is chosen in this 
particular combination. By definition this implies 



E(n,m) 
v) > =m 

i=i 

Moreover, we introduce the combinatorial sum symbol S-p^^^^ , 



(13) 



(14) 



m) 



where / is a function in the components of the vector i7("'™) , and where the sum contains 
all possible combinations with repetition for the situation that m elements are taken 
out of a set with n elements. E.g. it is obvious that for n = 2,m = 3 the equation 

f = /(3, 0) + /(2, 1) + /(I, 2) + /(O, 3) (15) 



■^^(2,3) 



is valid. Now the preparation for the virtual analysis is complete, so we can introduce 
the announced closed form of the wave functions \l/5(ri, k) 

k n—1 



(16) 



This is proven in Appendix A. The expression ( fT6ll allows an analysis of the norm given 
in Sec. HI 



Now we show that the simple algorithm is a one step method that applies pertur- 
bation theory correctly in first order in At for all /c G N. Therefore we expand the with 
perturbation theory calculated wave function ^(At, k) for all 6 N in second order in 
the time step At, where we regard as a starting point the Eqn. (|9]) for t = At: 

At 



^(At, A;) = e 



-iHoAt 



*(0) 



I I dt'e 





-^Hoi^t-t')W^t')^^t',k-l) 



Ho + WiO) 



At' 



+ HoWiO) + (1 - Ski) Wioy 



Hi + W{mo 



dW{t') 



dt' 



t'=0 



#(0) + C(At^) 



(17) 
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Then we calculate with Eqn. fill I) the wave function ^5(1,/;;) for all /c G N by 
propagation of the start wave function \I^(0) over one time step with the simple algorithm 
and expand the result in second order in At. For this calculation it is practicable to 
write W{At) instead of W{1): 

^s(l, k) = e-*^°'^*^(0) - iAtW{At)^sil, k-1) 
= e-*^°^*^(0) - iAtW{At)e-'^°^' x 

-rHoAt^^Q^ _ _ ^^^^ W{At)^s{l, k-2) 



X 



-iHoAt 



- iAtW{At)e 



-iHoAt 



-At^l-6ki)W{Atye 



-iHoAt 



+ C(At=^) 
1 - ^ ( ^0 



+ 2(1-5, 



1^(0) j At - 
W{0f + 2t 



Ae 
2 

dW{t' 



Hi + 2iy(0)^o 



dt' 



t'=0 



^ (0) + 0{At^ 



By comparison of (fT7|) and (fT8|) it can be cognized that 

At2 



*(At, k)-^s{l.k) 



dW{t') 



dt' 



t'=0 



^(0) + C(At=^ 



(19) 



so the simple algorithm is a one step method which applies perturbation theory correctly 
in first order in At for all A; G N. Due to standard numerical textbook analysis of the 
asymptotic development for the global discretisation error of one step methods [19], the 
difference between the with perturbation theory calculated wavefunction \l/(t, A;) and 
the wave function propagated over the time t with the simple algorithm \l/5(n, A;), with 
n = t/ At, is given b}l§|: 

^(t, A;) - ^shh k) = At x(t, A;) + ^(At^) (20) 

For the function x(t, k) appearing in the above equation holds that it is independent of 
the time step At and it fulfils x(0, A;) = 0. 

For the evaluation of the accuracy of the simple algorithm has to be taken into 
account that the with perturbation theory calculated wave function \I^(t, k) is itself an 
approximative solution of the Schrodinger equation ([2]), which deviates from the exact 
solution \l/(t). This deviation can be noted as \l/(t) — \l/(t, A;) = (f){t,k). Therefore for 
the difference between \l/s'(n, A;) and the exact solution \l/(t) holds: 

^ (t) - ^s{n, k) = 0(t, A:) + At x{t, k) + ^(At^) (21) 



§ The associated proposition in ^9] to the Eqn. ([20| was discussed there for real functions instead of 
a complex wave function but it is straightforward to see that this is no limitation for an application 
of this proposition here. Futhcrmore we suppose that the electric field E{t) and therefore the wave 
function \E'(t), too, depends in practical applications of the simple algorithm smoothly on time so that 
the proposition requirements given in |19j are not violated. 
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As the main conclusion of this section, the simple algorithm is a one step method that 
applies perturbation theory correctly for all /c G N only in first order in At, so one might 
think that it is not useful to calculate wave functions \E'5(n, k) with k > 1. However, this 
reasoning is not correct because the use of a higher perturbation order k takes according 
to ( I2T]) influence on the difference \E'(t) — ^^5(71., k) and stabilizes so the simple algorithm 
against divergences of the norm, see the appendant discussions in Sec. [51 

4. Error analysis of the norms for the wave function k) 

4-1. Preliminary remarks 

From the closed form expression ( fT6l) for the wave function \E's(?t,, k), it emerges that the 
wave function can be decomposed into terms \E'm ^(n) for different orders of perturbation: 

k 

^s{n,k) = J2^rn,s{^) (22) 

m=0 

with 

^m,s{n) = {-^^tr n(^(^-^T'"-'e-^^«^*)^(0) (23) 

where the parameter m denotes the order of perturbation. The wave functions "^m.sin) 
have a clear interpretation as they are related to an interaction between the laser pulse 
and the molecular system, where m photons are exchanged, and therefore the electronic 
state changes m times during the time interval [0,n At]. Employing the decomposition 
(|22l) . the norm yields: 

Ks= {^sin,k)\^s{n,k)) 

k k 

= Y.Y.{^^'S{n)\^HAn)) (24) 

j=0 h=0 

Using the substitution p = j + h, we can transform (124]) regarding that the norm of the 
initial wave function ^(0) is defined as ^^(0)|^(0)^ = 1 as follows: 

k j-\-k 

<5= EE(*:'-^(^)l^^-^-^(^)) 

j=0 p=j 

2k min(p,fc) 

p=0 j=max(0,p— fc) 
2k mm(p,fc) 

= 1 + E E (*.,5(n)|#p_,,5(n)) (25) 

P=l j=max(0,p— fc) 

We note that terms for odd p in (l25l) are zero, which results from the choice of the 
initial conditions, where only the electronic state |1) is populated (see (EI)). These terms 
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involve the scalar product of wave functions in the different electronic states |0), |1), 
which are orthogonal. Due to this connection we substitute p = 2m in fl25|) . which 
yields: 

k min(2m,fc) 

m=l j=max(0,2m-fe) 
k 

:= 1 + 5^ (26) 

m=l 

where the terms N^2mS^ given by 

m\n{2m,k) 

K2m,s= E {^jAn)\^2m-,An)), (27) 

j=max(0,2r?i— fc) 

characterize norm deviations from unity. As was already discussed in [1] and [15] these 
orders A^n,2m,5 can be decomposed into two different types: For 2k > 2m > k we call 
these orders oscillatory orders: 

k 

K2m,s= E (*.,5(n)|^2^_,,5(n)), (28) 

j=2m—k 

and for k > 2m > they are called stationary: 

2m 

K,2m,s = E {^UnW^m-Un)) . (29) 

Note that the stationary orders contain no explicit dependence on k any more. 
In [Ij, it has been shown that these two kinds of orders behave differently when 
parameters in the numerical simulation are changed and thus they have different physical 
interpretations. The former presented discussion was based on equations, which were 
presented without proof. In what follows the missing derivation will be given for the 
simple algorithm. 

4-2. Norm analysis for the simple algorithm S 

In the norm analysis for the simple algorithm in this chapter, we state first that the 
stationary orders can be written in a closed form. Secondly we discuss that for the 
oscillatory norm orders we can introduce an approximation that allows to analyse how 
these norm orders scale in the time step At, and in this context we explain an evidence 
called annihilation thesis, which is a pre-condition for the approximation. Thirdly we 
show that we can easily apply this approximation method, implemented foremost for 
the oscillatory orders, for the stationary orders, too. 

In the calculations given here we are focussing on the considerations referring to the 
mentioned approximation because the most important point for the understanding of 
the mathematical background of the results presented in [1] is to get the idea how this 
approximation is introduced. Other, more straightforward calculations are given in the 



Divergence of time- dependent perturbation theory 



10 



appendices. 



Thus, for the stationary norm orders A^^ 2m s perturbative wave functions cal- 

culated via the simple algorithm we first note that 



(n,m) 



(30) 



the proof of this equation is given in Appendix B We emphasize that in (!30!) only even 
orders of W{n) appear, for which ([3]) is valid. 



After having presented the closed form fl5Ul) for the stationary orders of the simple 
algorithm A^, 



n,2m,S'> 



k > 2m > 0, we will investigate now the announced properties of the 
oscillatory orders A^^ 2m,5; ^ < 2m < 2k. To do so, we employ ( |28!) to calculate that for 
the orders N^ 2ms^ ^ < < 2k the following equation holds: 



E 

j=2m.~k 



2m-j 



E E ' 

j=2rn-k j) 7'^(„,2m-j) 



*(0) 



n 

q=l 



W{q) 



X 



n-1 



[n 



{n,2m—j) 

pY"-p e~ 



^(0) 



and P= 



(31) 

2m-j) is a product of the 
The same bracket terms 



X n {^(^ 

Each summand of the sums in fl3T|) over j, V^^nj) 
sign factor (— 1)™"-' and a bracket term ^^(0) |- ■ \E'(0) 

^(0) | - ■ "I ^ (0)^ appear multiple in different summands of these sums, but the j-value 

of these summands differs, and thus they have different sign factors (—1 )'""■'. Due to 
these different sign factors, these summands must cancel (at least) partly each other, 
and for each of the various types of bracket terms in the sums in ( l3Ti) can only survive 
bracket terms with either the sign factor (—1) or (+1). Now it can be shown that 
all these surving bracket terms do not have different signs but the same sign, scilicet 
(— l)'^""^. This coherence is from now on called the annihilation thesis and the proof of 



it is presented in [Appendix C 



Our aim is to count the total number 2m s surviving bracket terms in (l3Ti) 

because the scaling of the number of these terms in the time step At helps to draw 
conclusions how the oscillatory order N^2m s scales in the time step At. Employing the 
annihilation thesis allows to count the total number of all the surviving bracket terms, 
and it can be concluded that: 



k 

71, 2m, 5 



E 

j=2m—k 



E^ E 



1 



(32) 
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Now we introduce an approximation for the calculation of the norm orders A^^ 2m s 

'. 2m 

approximating each of the surviving summands in (13T|) with a constant factor W 
and we take into account the fact that all the surviving terms have the sign-prefactor 
(— 1)'^"'" by a global sign prefactor (—1)'^"'" for the calculation of the norm orders 



j=2m~k 




(33) 



Since the sum is executed over all combinations with repetition for a choice of 



j elements out of a set with n elements, elementary statistics yield [20 



(n + j - l)(n + j - 2) ■ ■ ■ (n + l)n 




(34) 



Inserting ( l34l) in ( l33l) and shifting the sum index j by m, we find 

k—m 



n,2m.S 



._^Nfc-m^^2mp^2m^2m 



(-1)-' 



j=-{k-m) 



[m + i)\{m- i] 



0|i 



^At^^iy |^(m, A;) I 



where we introduced the help function ^(m, k): 

k—m 



{m + i)\{m — j)\ 



(35) 



(36) 



j=-(fc-m) 

With an induction proof over k one can show that if the oscillatory condition k < 2m < 
2k is valid, ^ (m, k) can be written in the form: 

-^\k—in 



^(m, k) 



(37) 



m k\{2m — 1 — k)\ 

As the base case for this induction proof we choose k = m, because for the validity of 
the oscillatory condition k < 2m < 2k this is the lowest possible value for k, and we 
calculate with (l36l) that C,{m,m) = l/(m!)^, which fits with the result we get for k = m 
with ( 137|) . Then we perform the induction step assuming as induction hypothesis that 
for an arbitrarily chosen value for k satisfying the oscillatory condition (1371) is true: 

k+l—m 

i{m,k + l)= Yl 



j=-{k+l~m) 



(m + j)l{m — j)\ 
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2 



(M~m 1 + 



{k ^ \)\(2m - k - \)\ {m + j)\{m - j)\ 

j — —{k—m) 



^' (fc + lj!(2m - - Ij! m fc!(2m - 1 - fcj! 

-| k+l—m -| 

(3J 



m + 1)! [2m- 1 - (A;+ 1)]! 

This is fl37l) for the shift k — )■ k + 1, so the induction step and the proof of (EZ]) is 
complete. 

For a re- writing of f l35p . we first substitute the expression (13 7p for ^(m, k) in fl35l) . and 
second, introducing the propagation time t = n At, we substitute furthermore for n 
in fl35|) . So as a final result we get Eqn. (14) of Ref. |T] for the simple algorithm: 

/I \k-m j-2m 

In the latter equation we have replaced the index n on the left side by t referring ex- 
plicitly to this time. 

As the third task in this chapter, the above approximation can be applied not only 
for the oscillatory orders N^2mS^^ ^ < 2k, but as well for the stationary orders 
^n,2m,sy ^ — > 0, because in (130|) it can be seen that all terms appearing in the 
sum over V^i^.m) have the same global sign (—1)™. Therefore we approximate all factors 
W{jf by PF^ in ([30]). Regarding dMD, this leads to 



= (-i)™At2-1F'" E 1 



(n,m) 
3 



O 



(40) 



Substituting n = we receive the final result 



Nt2m.,s = {-ir—At"^w'"' + O {At'--'') , (41) 
which is the Eqn. (13) in Ref. [1] for G = S. 

In Sec. [5] we will discuss the consequences of the results derived in this chapter. 



5. Interpretation 

In this chapter we will pick up the areas of interest established in Sec. 2.3 and Sec. 
3 of [1], explain the mathematical background to the discussions related to the simple 
algorithm we started there and carry them forward. In line with this task we give 
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Figure 1: Time-dependent norm of wave functions determined by perturbation theory 
(simple algorithm) of different orders k, as indicated. Curves obtained in different 
order k of perturbation theory are compared for time steps of At = 4 x 10~^ fs and 
At' = 2At = 8 X 10-2 fs. 



the results derived in Sec. S] a physical interpretation. Moreover, we present in a 
gedankenexperiment an analogy of perturbation theory to a scattering experiment. 

As a starting point we will discuss the stationary orders. These orders depend in 
leading order on the m-th power of the time step At. For a not too long propagation 
time t and a small time step At, the stationary order for m = 1 (calculated with ( HTl) ) 

iV/;2,5 = -t At + 0{At^) (42) 

gives the dominant contribution to the norm deviation. Thus, the norm deviations 
caused by the stationary orders are negative, and they are independent of the order k 
but they depend approximately linearly on the time step At. So, in the limit At — 0, the 
norm deviations caused by the stationary orders vanish. Thus the stationary orders are 
related to the At-dependent part in the difference of the approximative wave function 
"^si^jk) to the exact wave function \l/(t) (see Eqn. ( 12T1) ). From these findings we 
conclude that these errors are numerical errors which arise using the simple algorithm 
from the approximation of the integral in ( fTOl) by only one summand, and they have 
no physical meaning. This interrelation was discussed in the explanations to figure 3 in 
[1] , shown in this paper for clarity again in figure [H where one can realize first that for 
early propagation times a bisection of the time step At leads to a bisection of the norm 
deviations caused by the stationary orders, and second that these norm deviations are 
the same for equal time steps At but different orders k. 

According to ([3]) and ( 130|) the stationary orders N^^2m,s '^^^J depend on the time 
step At, the dipole moment /i and the electric field E{t), but they do not depend on 
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{k = 14) involved different potential gradients mo = Aj x 10~^ a.u., with numerical values 
of Ai = 3, A2 = 2 and A3 = 1. Whereas for the steepest potential the norm is conserved, 
divergences occur for smaller gradients where the wave-packet motion proceeds slower. 




the Hamilton operator Hq of the unperturbed system. This is one more evidence that 
the stationary orders are only numerical errors. Thus in particular they do not depend 
on the shape of the potential surfaces Vj{R),j G {0, 1} of the two electronic states |1) 
and |0), nor on the propagation of the wave function components \&i(-R, n) and '^o{R, n) 
(here we explicitly wrote out the -R-dependence) onto these surfaces. This context was 
presented in the discussion of figure 5 of [1], shown in this paper in figure [2] for the case 
when the simple algorithm is used: For a simulation model with linear potentials 

Vj{R) = nijE + Cj, J G {0, 1}, mi = -mo (43) 

we varied in |T] the parameter for the steepness of the potentials mo. As a result of the 
simulations it is obvious that for small propagation times the negative norm deviations 
caused by the stationary orders are the same. 

More explicitly, we compute the norm deviation caused by the stationary orders 
N^ 2m 5' ^ — > for an electric field given by 

E(t) = A(t)cos($(t)) , (44) 

where A{t) is the envelope of the electric field and <l>(t) is a fast oscillating phase 
function. Since we can assume that the lowest stationary order A'^^ 2 s gives the 
dominant contribution to the norm deviations, it is in our interest to calculate this 
order analytically. Therefore we consider that permutation vectors of the form z/*^"'^) 
have n components at which one component is 1 and all other components are 0. So we 
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can simplify f l3Up for m = 1: 



= - E n ^(^■) 

n 

= -At At (45) 

Inserting (|3]) and (jH]) into f l45|) and approximating the sum by an integral leads to the 
result: 

^i,5 = At /"* rft' A(t')' cos' ($(t')) (46) 

Here, we denoted the time dependency again by t. Since this integral in the limits and 
t is proportional to the stored energy in the laser pulse in this time interval, the norm 
deviations caused by the stationary orders at a certain point in time t are approximately 
a proportion for this quantity. Therefore for the limit t — j- oo, the stationary order N^2 s 
is proportional to the total energy of the laser pulse. Within the " slow- varying envelope 
approximation" [211 122] we can approximate ( 146|) by substituting the factor cos' {^{t')) 
under the integral by | and get as a result: 

K2,s = I dt' (47) 

This result reveals that if the "slow- varying envelope approximation" is valid, the norm 
deviations caused by the stationary orders do not depend approximately on the phase 
$(t) but on the time integral over the squared envelope of the electric field. Futhermore 
it can be cognized from (147]) that N^2 s depends linearly on the time step At. This result 
is related to the fact that according to Eqn. ( !2T1) the leading order of the At-dependent 
part of the deviations of the wavefunctions \l/5(n, /c) to the exact solution \E'(t) is the 
first order. 

In particular we regard as an example for the use of ( H7I) the Gaussian laser pulse 
modified by a linear spectral chirp 62 we employed for our numerical application of the 
simple algorithm in [1] and compare the results we got there numerically with the norm 
deviations we are now able to calculate analytically with (HTj) : 
The unchirped pulse is given by 

E{t) = E'^e-^'^'-''^^" cos [uo (t - t,)] (48) 



with an envelope having a full width at half maximum of t' = v/41n2//3', and the 
chirped laser pulse is given by 

E{t) = Eoe-^(*"*^)' cos [uo (t - t,) + I (t - t,)'] (49) 

In the last equations, the field strengths are denoted as Eq and Eq for the unchirped and 
for the chirped fields, respectively, and td denotes the point in time when the envelope of 
the field is maximal. The various parameters appearing in the equations for the shaped 
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Figure 3: The norm deviations arising from a numerical simulation employing the simple 
algorithm (panel (a)), and the norm deviations calculated analytically with f lSTjl (panel 
(b)). The analytically calculated deviations in panel (b) accurately match the numerical 
results in panel (a). 
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and for the unshaped electric fields are related as follows |23j : 

Eq 



E' 



a2 



/lfl/2 ' "9 



(50) 



4/3'2 I ^2 

For the electric field given by (149!) we can calculate with (147|) and (!50l) that the norm 
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deviations caused by the stationary orders are given by 



t,2,S 



TX 



128 In 2 



1 + erf 



32 In 2 r 



12 



r'4 + (161n2 62) 



{t - td) 



(51) 



Norm deviations calculated with (|5T1) match excellently with the numerical results which 
are presented in figure 4 of [1] . This relation is presented in figure [3l where the norm 
deviations from unity caused in the numerical simulation with the simple algorithm are 
shown again in panel (a), and the norm deviations from unity calculated analytically are 
shown in panel (b). For the analytical calculation we used the parameter values which 
lead to the numerical results presented in figure 4 of [1], namely fi = 1 a.u., At = 3.31 
a.u., t' = 4.13 * 10^ a.u. and -Eq = 1.19 * 10~^ a.u. and different values for the spectral 
chirp parameter 62 , which are given in figure [3l 

Furthermore, in the limit t — ?■ 00, one obtains from 



hm N, 



t—>oo 



t,2,S 



At E'i^ r' 



, • (52) 
321n2 ^ ' 

This result is independent of the spectral chirp parameter 62 because the total energy of 

the laser pulse does not depend on the spectral chirp parameter 62- For the parameters 

of our numerical example, (l52ll leads to a value of 0.927, which is in excellent agreement 

to our numerical results. 



Next, we discuss the norm deviations caused by the oscillatory orders A^^ 2m s^^ ^ < 
2k. For longer propagation times t these orders make up the dominant contribution to 
the norm deviations for the simple algorithm because in taking the ratio between the 
stationary and oscillatory orders (see (14T1) and (!39|l ). this number scales with (At/t)™ 
thus for large t the stationary terms are negligible. 

As was discussed in [1], the oscillatory terms, in leading order, do not depend on 
the time step At. Therefore we can conclude that the oscillatory orders N^2ms^^ < 
2m < 2k in contrast to the stationary orders N^2mS^^ — > do depend on the 
At-independent part of the difference of the approximative wavefunction ^5(77,,/;;) to 
the exact wavefunction \I^(t), which is given in Eqn. (^T^ by (f){t,k). Futhermore, for 
the At-independent part of the oscillatory terms on the norm deviations, the following 
argumentation holds: 

The factor m [k\ {2m — 1 — k)\] in the denominator of the ratio in (l39l) increases for an 
increment of the parameter m of an oscillatory order N^2mSi ^ ^ < 2k. However, 
the term t^™ in the nominator of these equations overcompensates this factor for large 
values of the propagation time t. Thus, the larger the order parameter m is, the larger 
has to be the propagation time t, in order that the contribution of the corresponding 
oscillatory order N^2mSi ^ ^ 2m < 2k leads of all oscillatory orders to the largest 
modification of the norm deviations form unity. Due to the factor (— 1)^~"* in f[39|) . this 
modification of the total norm deviation is negative or positive. Moreover, the sign of 
the terms alternates as a function of m, so that with increasing propagation time t, the 
sign of the contribution of the oscillatory orders to the norm deviations changes. This 
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is the actual reason why we named these terms oscillatory orders in [T]. 

Eventually, for large propagation times t, we get due to the positive sign of the highest 

oscillatory order (m = k) 

,2k 

N^s = JT^W"' + O (At) (53) 
(A;!) 

a divergence towards +oo, since the exponent 2m in ( 139|) reaches its biggest value for 
m = k. Nonetheless, we can retard or even suppress the point in time when the norm 
deviations from unity caused by the oscillatory orders become a relevant contribution 
by using a higher value for k in the simulation. This retardation can be seen in figure [H 
where the use of /c = 14 instead of /c = 6 leads to a later point in time for the divergence. 
This effect can be interpreted as follows: 

The order k represents the maximal number of photons, which interact with the molecule 
in the numerical simulation. Since the norm deviations caused by the oscillatory terms 
in leading order do not dependent on the time step At but on the order k (see ( 139|) ). 



the norm deviations from unity caused by the oscillatory orders N^2mSi ^ ^ 2m < 2k 
are related to the truncation of the expansion of the wave function ^^5(77., k) for higher 

th 

orders than the -order of the multi-photon processes. Since for the simple algorithm 
both norm deviations related to the stationary and to the oscillatory orders appear, 
we have to differentiate in the use of the simple algorithm between these two types of 
contributions to the norm deviations. 

As a next task, we show that from the results presented in Sec. |4]we can conclude 
that the oscillatory orders N^2m si'^^ — ^ ^ contrast to the stationary orders) 
depend on the shape of the potential surfaces Vj{R). In order to explain this, we 
calculate explicitly with ( l23l) and ( l28l) the oscillatory orders for the simple algorithm: 

k I 

K2m,s= (-i)™-^At2'" 5^ J2 (^(0) 

n-1 



I X 



9=1 



#(0) , (54) 



p=0 

The bracket terms (\E'(0) | - ■ - j \E'(0)) in (JSlj) describe the overlap between a bra-state that 
is influenced by interaction operators W{q) at certain points in time ^i^^"'''^ 7^ with 
a ket-state that is influenced by interaction operators W{n — q) at other points in time 
^pi^Iq™""''' 7^ j . Through the impact of the interaction operators at different points in 
time the bra-state and the ket-state propagate over the time interval [to, t„] in a different 
way. 

That means for points in time t„/ with to < t„/ < t„ the bra-state and the ket-state 
can be localized in different electronic states. But at t„, when the overlap is calculated, 
the bra-state and the ket-state must be in the same electronic state, otherwise the 
overlap is zero. If the absolute value of the difference between the potential gradients in 



the two electronic states AV'{R) :- 



dVi (R) _ dVo{R) 
dR dR 



is large in the spatial region, where 
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the wave function is positioned during the interaction of the molecular system with the 
laser pulse, the differences in the temporal propagation in these two electronic states 
|1) and |0) lead to a small overlap of the bra-state and the ket-state. This leads to the 
effect that the oscillatory orders are suppressed for large potential gradient differences 
AV'{R). For this effect we presented already in [1] numerical evidence without analytical 
proof. In order to comprehend this evidence, we look at figure [2] again: For an increase 
of Ay (i?) = 2mo (see f H3|) ) the divergence related to the oscillatory orders is retarded 
temporally, and for the largest value mo = 3 * 10~^a.u. the divergence can even be 
suppressed. This effect can interpreted in that way: For a large potential gradient the 
interaction time for the wave packet in the electronic state |1) with the laser pulse for a 
transition to the electronic state |0) is small and then this interaction can be discribed 
by a smaller amount of photons. As a result, in contrast to the stationary orders not 
only the electric field E{t) but as well the full dynamics of the molecular system are 
relevant for the values of the oscillatory orders N^2m ^A; > 2m > k. 

Having discussed the dependence of the oscillatory orders on the potential surfaces 
Vj{R), j G {0,1}, the dependence on the electric field E{t) leads to interesting phe- 
nomena, as well. Here in this context, we mention an effect that we discussed in [1] in 
detail: 

If the electric field E{t) is modified under the condition that the total energy of the 
laser pulse remains unchanged, what leaves the norm deviations caused by the station- 
ary orders unchanged in the limit t — )■ oo in good approximation, this is not true for 
the oscillatory orders. Namely, we showed in [1] that the energy conserving temporal 
broadening of the laser pulse by the introduction of a spectral chirp 62 leads to an en- 
largement of the norm deviations caused by the oscillatory orders and advantages norm 
divergences. This means that one needs because of this broadening of the Gaussian 
laser pulse processes where more photons are exchanged between the laser pulse and 
the molecular system in order to describe the interaction of the laser pulse with the 
molecular system in an adequate accuracy. Thus, this example reveals how the results 
explained here and in [1] can be used to gain deeper insight about the interaction of a 
molecular system and a laser pulse. 

Let us complete the discussion above by the following consideration: Here, our aim 
is to illustrate why the truncation of the perturbative expansion, taking only k orders 
in (I22I) into account, is a valid ansatz for an approximation, although this expansion de- 
composes the wave function into terms, which are not orthogonal (^(^^p,s\^q,s^ ^pq^- 
This illustration is done by a comparison of our situation with a gedankenexperiment, 
where a wave function is also split into non-orthogonal components. 
Therefore we imagine a transmission grating, that is irradiated with an electron beam, 
which has a Gaussian density profile. Due to the transmission grating, the wave function 
of the electrons after the transition through the grating is split into different compo- 
nents: Namely each of the slits creats a circular wave function, and the amplitude of 
each circular wave is related to the position of the according sht in the Gaussian density 
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profile before tlie grating. Tlie different circular wave functions interfere, and befiind 
tfie grating is a screen where tlie interference pattern caused by the electrons hitting 
the screen is monitored. Then we arrange an aperture in front of the transmission grat- 
ing, which is centered at the density profile of the electron beam. Now we adjust the 
minimal opening of the aperture, so that due to our accuracy of measurement we detect 
no difference of the interference pattern at the screen between the situation with the 
aperture and without it. So the result can be interpreted as follows: 
The interference pattern is created in good approximation both for the presence and 
the absence of the aperture by electrons which pass simultaneously through all the slits, 
which are for the presence of the aperture not masked by it, and we can conclude that 
because of the fast convergence of the Gaussian profile to zero in good approximation 
the electrons never pass all the other slits. So it makes no difference for the inter- 
ference pattern if we mask them by the aperture or not. For the construction of the 
pattern at the screen, we must not regard the different circular waves of the different 
slits independently of each other, but we have to take interference effects into account. 

In the same sense we can say that if we need k orders to suppress the norm deviations 
related to the oscillatory orders in our calculations, in the interaction of the laser pulse 
with the molecular system 0, 1, 2, . . . , photon processes happen simultaneously but in 
good approximation processes which need more than k photons do not happen. For 
such an interpretation of the norm deviations related to the oscillatory orders, we have 
to take into account that the different perturbation orders of the wave function are not 
orthogonal and interfere like the circular wave functions of the slits in the scattering 
experiment; in particular we cannot give an assertion with what probability V{x) a 
specific amount x of photons with < a; < A; participates in the laser pulse-molecule 
interaction. 

The norm deviations related to the stationary orders are not relevant for this 
consideration because they are only norm errors related to the numerical inevitable 
discretization of time. 

6. Summary 

In [1] we presented methods for the numerical application of perturbation theory with 
the aim to characterize the interaction of shaped laser-pulses with molecules. In order 
to control the quality of the calculations we analysed norm deviations caused by the 
necessary discretisation of the appearing time-integrals and also the truncation of the 
perturbation expansion; both causes can lead to substantial deviations. Moreover, in [1] 
numerical results on a model system incorporating a single nuclear degree of freedom 
were presented, where a chirped shaped laser pulse induces electronic transitions. 
These results, in particular the dependence of the different contributions to the norm 
deviation on parameters like the propagation time step, the steepness of the potential 
curves and the chirp parameter are explained in [1] for brevity by some equations and 
interrelationships without proofs. The latter are provided in the present paper for what 
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we called simple algorithm in [T], where the time- integral over one time step occurring 
in perturbation theory is approximated by one term. 

In our analysis of norm deviations from unity calculated with the simple algorithm in 
[1], we stated that two classes of terms of different character contribute to the norm 
deviation: The first ones are called the stationary orders and are of purely numerical 
nature. They can be suppressed in the limit of small time steps. The second kind of 
contributions, called the oscillatory orders, are related to the property of time-dependent 
perturbation theory, which is not norm conserving. These orders can cause oscillations 
in the norm of the total wave function, and moreover, for long enough propagation times, 
these terms can lead to divergences of the norm towards infinity. In this publication, 
we proof equations and interrelations which explain this behavior of the stationary 
and the oscillatory orders for the simple algorithm that was documented for numerical 
calculated results in [1]. Moreover, here we present a method to calculate the norm 
deviations caused by the stationary orders quantitatively analytically. The accuracy of 
this method is demonstrated for a numerical example. 

In [1] we stated furthermore that the oscillatory terms directly correlate with the 
order of the multi-photon transitions to be described, and by increasing the order of 
the perturbation theory the norm deviations can be reduced. In this publication we 
generate, beyond that, an easy interpretable graphic image of this reduction. Therefore, 
we compare our situation to the situation that a transmission grating is irradiated by 
an electron beam with a Gaussian density profile, what causes an interference pattern 
on a screen behind the grating. Moreover, there is in front of the grating an aperture 
centered around the electron beam, that we open further until the interference pattern 
does not change anymore. 

In the future, we will publish in |T3] how the calculation done in this paper can 
be devolved from the simple algorithm to what in [1] we called improved algorithm. 
Moreover, our research goals are to implement the results depicted in [1] and in this 
publication in order to get a better understanding of multi-photon processes. With 
the help of the derived results we are now in the position to analyse which processes 
of different orders are relevant for the correct description of a chemical reaction. 
Therefore we will compare the results calculated with converged perturbation theory 
and numerically exact results, which contain all perturbation orders. 
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In this appendix we prove that the wave functions \I/5(n, k) for the simple algorithm 
can be written in the form of ([T6l) . which is 



n-l 



{n,m) 



where we use the definition ( II ip for the simple algorithm as a starting point for the proof. 



As an inception of the proof we calculate as helpful lemmas with (fTT|) the following 
equations for "^si^^k) and \E'5(n,0): 

^5(1, k) = e-*^°^*^(0) - I At W{l)^s{l, k-1) 



m=0 
k 







m=0 



^(0) 



(Al) 



and 



n-l 



{n,m) 



(A2) 



For the proof of f|T6|) . we work with complete induction and the induction hypothesis 

q-p 



^s(p + 2,g-p- 1^ 



-iHoAt 



jn=0 'Pp(p+i,m) i=0 

for ag e Nandp = 0,1,2, ...,g- 1. (A.3) 

For p = and an arbitrary g G N is the induction hypothesis flA.Sp true because of 
f lA.ip . This implies that (lA.SP is true for the base case q = I- 

The induction hypothesis ( \A.3\\ implies by substituting phj p + 1 the equation 

'e h ^tr E n (^(p + 2 - jr^+™e- 

m=0 'Pp{p+2,m) i = 

for ag G N andp = 0, 1,2, - 2. (A.4) 

Because of (lA.2p we can assume that (lA.4p is correct for p = g — 1, too. With the 
definition f lTT]) for the simple algorithm we get 

^s{p + ^,q-p) = e-*^o^*^(p + 1, g - p) - i At (p + 2)^sip + 2, g - p - 1) 

for ag G Nandp = 0,1,2, ...,g- 1. (A.5) 



#(0) 
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iHoAt 



m=0 
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J2 H At)"+^iy(p + 2) ^ Y[ (W{p + 2 - j)"X2-r'e- 

^p{p+2,77i) i=o 



iHoAt 



m=0 



m=0 



,(p+2,m-l) 
P+2-i ( 



-iHoAt 



M/(0) 



(A.6) 
(A.7) 



g-p p+l 

^ (-Z At)™ ^ n (^^^ + 2 ~ jTp'+'^-r' e-*/^o At ' 

^m=0 'Pp{p+2,m) i=0 

for a g G N and p = 0, 1, 2, g — 1. 

In this calculation we regarded the fact that in flA.6p the term in the upper line contains 
only terms without the operator W{p + 2) and the term in the lower line contains only 
terms where W{p + 2) appears at least in first power in order to see that the terms in 
(1A.6P and ( 1A.7P are identical. Now we substitute in (IA.7P p by p — 1 and get: 



<?+i-p 



for a g G N andp = 1,2, (A.8) 

Moreover (lA.ip implies that f lA.SP is true for p = 0, too and this means that we have 
made the induction step and f lA.Sp is true for all g G N. With the substitutions n = p + l 
and k = q — p in (1A.3P the formula for the wave functions \l/5(n, fc) for the simple 
algorithm ( !T6|) is proven. □ 
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Appendix B. Proof of ( I30p for the stationary orders of the simple algorithm 

In this appendix we prove that for the stationary orders of the simple algorithm fl5U]) 
holds: 



{n.m) 



For the proof of f l30p we work with complete induction over n for an arbitrary m that 
fulfils the condition k > 2m > 0. 



As a proof for the base case n = 1, we derive with ( 129|) that for the norm orders 
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(B.l) 



Then as an induction hypothesis we suppose that all norm orders N^2m 5 for < r/ < 1 
with an arbitrarily chosen n G N suffice ( l30|) and show as induction step that this 
implies that A'!^+i^2m,s suffices ( 130|) . too. For this purpose we derive with ( p9|) ffist that 
for ^n+i,2m,s holds 

2m 

^n+l,2m,S = (^^-.^(^ + l)l^2m-i,5(ri + 1)' 
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Now we write out the operator W{n + 1) explicitly 

2m j 2m— j 

^^n+l,2m,S 



j=0 g=0 f=0 Pp(„j-9) 7'^„,2m-j-/) 
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9=1 



X 
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and then we introduce the new sum indices s '■= g + f and d := j — g: 

2m j 2m— d 
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Now we think about (lB.4p and cognize that there appears the threefold sum 
^j=o X]d=o Yl^I2j-d- -^'-'^ ^ threefold sum of this form and an arbitrarily function d, s) 
the equation 

2m— d 2m. 2m— s d+s 



2m 



j=0 d=0 s=j-d s=0 d=0 j=d 



(B.5) 



holds. The Eqn. (1B.5|) can be visualized by the picture that all 3-tuples {j, d, s) 
that are inside or on the boundary of a pyramid which is given by the four points 
(0, 0, 0), (0, 0, 2m), (2m, 0, 2m), (2m, 2m, 0), represent exactly all the combinations for 



j, d and s that are contained both in the threefold sum ^^™o Sd=o J2t=j-d ^^^^ 



\2m—d 



side of (1B.5P and in the threefold sum Yl^I^o Sd=o Sj=d right side of ( IB.5P 



thus, these sums are equal. The pyramid is presented in figure 
visualization can explained in detail as follows: 



for 2m = 8. The 



For the visualization of the threefold sum on the left side of ( IB. 51) it is helpful to 




Figure Bl: The pyramid, that we use for the visualization of the threefold sums 
^■=0 Ld=o 2^s=j-d on the left side and accordmgly l^^^^ l^.^^ on the right 

side of (lB.5p . is presented here for 2m = 8. Visualization means as explained in 
the text that all 3-tuples (j, d^ s) that are inside or on the boundary of this pyramid, 
which is given by the four points (0,0,0), (0,0,8), (8,0,8), (8,8,0), represent exactly 
all the combinations for j, d and s that are contained in these both threefold sums 
Ei=o Ed=o Es=i-d and Y.l=o Ed=o Yfjtd- Moreover, we have drawn in this figure the 
projections of this pyramid on the (j, (i)-plane and the (s, (i)-plane, which are triangles. 
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Figure B2: The triangle, that we use for the visuahzation of the outer twofold sum 
Sj=o ^d=o 3 ^^"^ ^^^^ ^^"-^^ fIB.Sp . is presented here for 2m = 8. We 

obtain this triangle by a projection of the pyramid in figure IBll on the (j, (i)-plane. 

think first about the projection of the pyramid on the j — plane; this is a triangle, 
which is given in figure [B2] for 2m = 8: 

All 2-tuples (j, d) that are inside or on the boundary of this triangle represent 
exactly all the combinations for j and d that are contained in the twofold sum 
Yl'j=o Yld=o- When we combine this result with the question which values for s for 
a given 2-tuple (j, d) inside or on the boundary of the triangle shown in figure [B2] induce 
3-tuples (j, d, s) that correspond to points located in or on the boundary of the pyramid 
in figure [BT| we can comprehend that all 3-tuples {j,d,s), that are inside or on the 
boundary of this pyramid represent exactly all the combinations for j, d and s that are 
contained in the threefold sum Yl'j=o J2d=o Yl1=j'-d ^^^^ ^^^^ °f ( 1B.5I) . 

By an analogue approach for the visualization of the threefold sum on the right side 
of ( IB. 51) it is helpful to think first about the projection of the pyramid on the s—rf— plane; 
this is a triangle, which is given in figure [B3] for 2m = 8: All 2-tuples {s,d) that are 
inside or on the boundary of this triangle represent exactly all the combinations for s and 
d that are contained in the twofold sum Yl1=o Yl'd^o^- When we combine this result with 
the question, which values for j for a given 2-tuple (s, d) inside or on the boundary of the 
triangle shown in figure IB31 induce 3-tuples (j, d, s) that correspond to points located 
in or on the boundary of the pyramid in figure IBlt we can comprehend that all 3-tuples 
{j,d,s), that are inside or on the boundary of this pyramid represent exactly all the 
combinations for j, d and s that are contained in the threefold sum Yl1=o Yl^d^o^ Yl'j=d 
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Figure B3: The triangle, that we use for the visuahzation of the outer twofold sum 
Ss=o Y^d^o^ o'^^^ -5 ^iid d on the right side of flB.51) . is presented here for 2m = 8. We 
obtain this triangle by a projection of the pyramid in figure IBll on the (s, (i)-plane. 



on the right side of (IB.Sp . 

Now for ( IB. 41) the summands ({j,d,s) in the threefold sum over j,d and s are of the 
special form 

aj,d,s) = {-iy4>{d,s) (B.6) 

and moreover, it is to show straightforwardly that 
q 

= (-irVp)W2,o, (B.7) 



Q 



(q—p) mod 2,0- 



(B.8) 



In the following calculations in this paper we will use both (]B.7p and ( IB. 81) . Thus, for 
a function ({j, d, s) for that (]B.6|) is valid, one can conclude with (IB.SP and (IB.7|) that 
holds 

2m j 2m— d 2m 2m— s d+s 

EE E (-i)w,^) = E E E(-iM^'^) 

j=0 d=0 s=j-d s=0 d=0 j=d 

2m 2m— s 

= E E ("^^'^ "'"'^ 2,0 Hd, s) 

s=0 d=0 
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With (IB.gp we can simplify ( ]B.4p and get as an intermediate result: 

m 2m— 2r 

^^n+l,2m,S ~ 
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(B.9) 



X ( ^(0) 



'■—0 d=Q 'Pp(n,d) 'P^n,2m-d-2r) 



9=1 

n-1 



X 



■ r"r A , T — r / (n,2m — d — 2r) 

p=0 

Reminding that W{n + 1)^*^' commutates with all operators, we simplify ( IB. 101) : 



*(0) 



(B.IO) 



^^n+l,2m,S 



r=0 



2{m-r) 



X (^(0) 



n 

9=1 



X 



n-l 



X 



p=0 



lW{n-py"-p e 



-iHoAt 



^(0) 



n,2{m-r),S 



(B.ii; 



r=0 



Employing the induction hypothesis, we can substitute now for the terms 
(IB. lip the right side of fl30|l and implicate that f l30|) is true for n + 1: 



2{m-r),S 



m 



m—r) 



r=0 
X 



E n^(^') 

'^0{n,m — r) 3 ^ 



2 1/ 



(n,m — r) 



n+l 



\2 u. 



(n+ l,m) 



^^(n + l,m) J = l 

Thus, the complete induction proof is succeeded and fl30l) is proven. □ 



(B.12) 
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Appendix C. Proof of the annihilation hypothesis for the oscillatory orders 
of the simple algorithm 

In this appendix we prove the annihilation hypothesis for the oscillatory orders of the 
simple algorithm, N^ 2m 5; ^ < 2m < 2k. 



As a starting point for this proof we take (13T]) with the substitution n — )■ n + 1 



2m,S 



-1 



X 



j — ^m — kVpl^n + l.j) 'P^n+l,2m-j) 



X ( ^(0) 



n+l 



(n + 1.2m— ?) . r"r a _i 



-p e 



#(0) 



J] (^e^^°^*iy(g)'^^ j J] (iy(n + 1 - p)^ 

q=l p=0 

and transform it in the same manner we transformed in Appendix B (1B.2P into (IB.4p . 
Thus, we get as a result: 

k j 2m— d 



-D-At- E E E E E (-1)' X 

j = 2m.~k d=0 S=j~dV^(ri,d) 'P^n,2m-d^s) 



X ( ^(0) 



n 

9=1 



n-1 



X W{n + iye 



s -iHoAt 



p=0 



n/ ^ {n,2rn, — d — s) 

(W{n-p)P"-p e 



-^_ffoAt^ 



^(0) 



(C.l) 



Now we think about ( 1C.1I) and cognize that there appears the threefold sum 

Yl^=2m-k J2d=o Yll=i-d- -^o^ ^ threefold sum of this form and an arbitrary function 
d, s) the equation 

k j 2m~d 2m min(2m—s, k) min{d+s, k) 

j=2m—k d=0 s=j—d s=0 d=max{2m—k—s, 0) j=max{2m—k, d) 

holds, thus we can transform fIC.lll into 



(C.2) 



2m min(2m—s, k) min{d+s, k) 

s=0 d=max(2m—k—s, 0) j=max{2m—k, d) 



X E E {m 

'Pp(n,d) 'Pp{n,2m-d-a) 



n 

9=1 



n-1 



X 



p=0 



^(0) 



(C.3) 



The Eqn. flC.2l) can be visualized by the picture that all 3-tuples {j,d,s), which are 
inside or on the boundary of a frustum of the pyramid discussed above that is the 
result of a truncation of this pyramid by the two planes j = 2m — k and j = k, 
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Figure CI: The frustum of the pyramid presented in figure IBlt that we use for the 
visuahzation of the threefold sums Yl^j=2m-k X]d=o Y^^J-d ^^^^ ^^"^ accordingly 

y^?'"n yymin{d+s k) right slde of flC.2D . is shown in this figure 

for 2m = 8 and k = 6 (the edges of this frustum are the highlighted ones). Moreover 
we have drawn in this figure the projection of this body on the (s, (i)-plane. 



represent exactly all the combinations for j, d and s that are contained both in the 
threefold sum Yl'j=2m-k X^d=o Yl1=j~-d ^^^^ ^^^^ ( 1C.2I) and in the threefold sum 

E2m \-^min(2m—s, k) ir-^minid+s, k) ,1 • t j_ • j c Art j-T 
n > . /o / n\ > • ; j\ on the right side of U.2 . thus these sums are 
s=0 ^d=max(2m—k—s, 0) ^—^j=max{2m—k, d) ^ !" 



equal. The frustum of the pyramid is presented in figure ICll for 2m = 8, A; = 6. The 
visualization can explained in detail as follows: 

The visualization of the threefold sum on the left side of f lC.2p is trivial, because it 
is easy to see by comparison of the left sides of ( IB. 511 and (IC.2I) that the truncation 
of the pyramid by the planes j = 2m — k and j = k only leads to a change of the 
summation limits of the outer j-sum. The visualization of the right side of ( ]C.2p is 
more complicated: 

Therefore it is helpful to think first about the projection of the frustum of the pyramid 
on the s — d— plane; this projection is presented in figure [C2] and it implicates that all 2- 
tuples (s, d) that are inside or on the boundary of this projection represent exactly all the 
combinations for s and d that are contained in the twofold sum Yll=o Yld^ax(2m-k-s o)- 
When we combine this result on the one hand with the question which value has the 
lower summation limit jmin for the innermost sum over j on the right side of (1C.2P for 
a given 2-tuple (s, d) inside or on the boundary of the projection shown in figure IC21 
we realize because of the shape of the frustum of the pyramid shown in figure 
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Figure C2: The polygon, that we use for the visuahzation of the outer twofold sum 

Z^s^o Z]d^ma^2m-fc-s, 0) ^'^^^ ^"^^ ^ "S^t side of flC.2p is presented here for 

2m = 8, /c = 6. We obtain this polygon by a projection of the frustum of the pyramid 
shown in figure [CT] on the (s, (i)-plane. 



that d < 2m — k implicates jmm = 2m — k and d > 2m — k implicates jmin = d (this 
result is visualized in figure [C3|) . On the other hand, when we combine this result with 
the question which value has the higher summation limit jmax for a given 2-tuple (s, d) 
inside or on the boundary of the projection shown in figure [C2l we realize again because 
of the shape of the frustum of the pyramid shown in figure [CT] that d-\- s < k implicates 
jmax = d + s and d + s > k implicates jmax = k (this result is visualized in figure IC4I) . 
Thus the innermost sum over j on the left side of (]C.2p is of the form Yl^l^tx{2m-k d)- 
Now for (IC.ll) the summands C(j, d, s) in the threefold sum are again of the special 
form flB.6P and with (IB.7P and fIB.SP we get in a fourfold case differentiation the following 
results for the calculation of the innermost sum over j on the right side of (IC.2p . where 
we name the four cases for afterwards comprehensible systematical reasons An+i-, -Bn+i, 
Cn+i and Dn+i- 

An+i'- d < 2m — k and d + s < k: 

min{d+s, k) d+s 

J2 {-iy^{d,s)= Yl (-l)'<^(^'^) = (-l)'Wfc)mod2,o'^K^) (C-4) 

j=max(2m—k, d) j=2m—k 
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Figure C3: The value for the lower sum limit jmin for the innermost sum over j on the 
right side of ( ]C.2l) depends on the values for the summation indices s and d of the two 
outer sums of this term and moreover on the parameters 2m and k: For d < 2m — k 
holds jmin = 2m — k and for d > 2m — k holds jmin = d. The two different areas in the 
the (s, (i)-plane for the different kinds of dependencies of jmin on s and d can be derived 
by the form of the frustum of the pyramid shown in figure IClj they are illustrated in 
this figure for 2m = 8 and k = 6. 




Bn+i'- d < 2m — k and d + s > k: 



min(d+s, k) 

E (- 

j=max{2m'~k, d) 


-iy^{d,s) 


k 

i=2m—k 


= (-l)V(c?,s) 


(C.5) 


Cn+i'- d > 2m 


— k and d 


+ s > k: 






min(d+s, k) 

E (- 

j=niax{2m—k, d) 


-iy<p{d,s) 


k 

= 5^(-l)W,^) = 

j=d 


(-l)'Vrf) mod2,o<^(^'^) 


(C.6) 


Dn+1- d > 2m 


— k and d 


+ s < k: 






min(d+s, k) 




d+s 







j=max{2m—k, d) j=d 

In figure IC5I we marked the different areas in the s — d-plane for the four different 
cases An+i, -Bn+i, C'n+i and -D„+i and the parameter values 2m = 8 and k = 6. 
Because in (]C.3P appears an sign factor (—I)'", which is not included in the fourfold case 
differentiation, and in the three cases C„+i we get an additional sign factor 
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Figure C4: The value for the higher sum hmit jmax for the innermost sum over j on 
the right side of (IC.2|) depends on the values for the summation indices s and d of the 
two outer sums of this term and moreover on the parameters 2m and k: For d + s <k 
holds jmax = d + s and for d + s > k holds jmax = k. The two different areas in the 
(s, (i)-plane for the different kinds of dependencies of jmax on s and d can be derived by 
the form of the frustum of the pyramid in figure IClt they are illustrated in this figure 
for 2m = 8 and k = Q. 




(—1)'^, we have shown for the three cases -B„+i and C„+i that the corresponding 

terms contributing in the sums on the right side of fIC.SP are proportional to a global 
sign factor (—1)'^"™. 

Thus we have proven that these terms fulfill the annihilation thesis, so for us only 
remains the task to show that the terms contributing in the sums on the right side 
of (IC.3|) and corresponding to the case -Dn+i, d > 2m — k and d + s < k, fulfill the 
annihilation thesis, too. 

For the solution of this problem it is practicable to call the sum over all terms on 
the right side of (]C.3p . which fulfill the condition for the case -D„+i, in our following 
calculations -D^*^|2mS- this notation, the superscript index tn+i in the brackets 
denotes that the propagation time t is equal to (n + 1) At, while the subscript + 1 is 
related to the fact that the case differentiation above, in which the case -Dn+i appears, 
is required in the context of a calculation where we wrote out the operators W{n + 1) 
appearing in the bracket term in ( 13T|) explicitly (see the transformation form ( 13T|1 to 

(EH)). 

We can note these terms D^nVi^'hnS with the conditions d > 2m — k and d + s < k 
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Figure C5: The corresponding areas inside the polygon presented in figure IC2I for the 
four different cases -B„+i, C„+i and D^+i presented in the text are shown in 

this figure for 2m = 8 and k = Q. (The points being located on boundaries bordering 
on different areas can arbitrarily allocated to one of these areas allocated cases; the 
assignment chosen for the definitions of the four cases in the text was made because it 
makes the calculations easier.) 




and (1C.7P in the following way, where we regard that the inequations s < k — d < 
k — (2m — k) = 2{k — m) hold: 

min[2(k—m) ,2m\ min(2m—s, k, k—s) 

^?;t2tU = (-l)-At^™ E E (-1)'^. mod 2, X 



X 



X 



E E 

'Pff{n,d) 'P p{n,2m-d-s) 



s=0 d=max{2m—k, 2m— k—s, 0) 

n 



n-1 



q=l 



p=0 



^(0) 



(C.^ 



Now we consider that for the oscillatory orders holds k < 2m < 2k, thus we can simplify 
the max and min functions in the sum limits in the sum over s and d in ( IC.SP and taking 
(jS]) into account we get with s substituted by 2r as a result: 

k—m k—2r 



r=0 



d=2m-k Vp(n,d) 'P^n,2{m.-r)-d) 
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vl>(0)) (C.9) 



As a next step we rename first d by j and second we write the operator W{n) explicitly 
out in the same manner like we did this in (IC.ip and flB.4p for the operator W{n + 1): 



k—m 



D 



(t„+i,k) 
n+l,2m,S 



r=0 



'Pp(n-l,d) 'P^7i~l,2(m.-r)-d-s) 



k-2r j 2{m-r)-d 

E E E 

j=2m—k d=0 s=j—d 
n-1 



(n-l.d) 



X 



n-2 



X iy(n)'e 



p=0 



(C.IO) 



Then we think about (ICIOP and realize that there appears the threefold sum 
X]j=2m-fc X]d=o Xls^-d^"'^- -^o^ ^ threefold sum of this form and an arbitrary function 
d, s) the equation 

fc— 2r J 2(m.— r)— d 2{m—r) min[2{m—r)—s, k—2r] min{d+s, k—2r) 

c(j,rf,^)= EE E (c.ii) 

j=2m—k d=0 s=j—d s=0 d=max(2m—k—s, 0) j=max{2m—k, d) 

holds. We can prove (IC.llI) by substituting the variables k hj k — 2r and 2m by 
2(m — r) in flC.2p in the upper sum limits. The reader can see for oneself that one can 
also visualize both sides of ( IC.llI) in a similar way we did this for ( ]C.2I) before. With 
( IC.llI) we can transform ( 1C.10|) into 

k—m 

D(*::nl), = f-l)-At2-^ty(n + l)2- X 



'n+l,2m,S 



X 



r=0 



2{m—r) min[2(m—r) — s , k—2r] min{d+s, k—2r) 

EE E (-i)'x 

s=0 d=max(2in—k—s, 0) j=max{2m—k, d) 

n-1 



X 



E E 



*(0) 



n 

q=l 



n-2 



{n — l,2(m — r) — d — s) 
7i — l — p f 



^iHoAt ^ 



-iHoAt 



(C.12) 



X W{ny e-*^o^* Yl [W{n - 1 - pY 

p=0 

and again we get with (1B.7P and (IB.SP in a fourfold case differentiation the following 
results for the calculation of the sum over j on the right side of (IC.lip . where we now 
name the four cases An, Bn, Cn and 
An', d < 2m — k and d + s < k — 2r: 

min{d+s, k—2r) d+s 

E {-iy<P{d,s)= E (-l)'<^(^'^) = (-l)'W+fc)mod2,o<^(^'^) (C-13) 

j=max{2m—k, d) j=2m—k 
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Bn'. d < 2m — k and d + s > k — 2r: 

min{d+s, k—2r) k—2r 

(-WK^)= E (-i)W^) = (-i)VK^) (C.14) 

j=max{2m—k, d) j=2m—k 

Cn'- d > 2m — k and d + s > k — 2r: 

min{d+s, k—2r) k—2r 

(-l)W^) = Yl(-iy^id,s) = (-l)V'i)mod2,o<^K^) (C-15) 

j=max{2m—k, d) j=d 

Dn'- d > 2m — k and d + s < k — 2r: 

min{d+s, k—2r) d+s 
j=max{2m—k, d) j=d. 

In an analogue argumentation as in the discussion of (IC.Sp . in (]C.12p a sign factor 



(— 1)*" appears, which is not included in the calculations done in the fourfold case 
differentiation. Since in the three cases y4„, C„ we get an additional sign factor 
(— 1)'^ like for the three cases An+i, Bn+i, Cn+i in the discussion above, we have shown 
for the three cases 5„, C„ that the corresponding terms contributing in the sums 
on the right side of ( 1C.12P are proportional to a global sign factor (—1)^""*. 



So we have proven that these terms fulfill the annihilation thesis, and now, the task 
that remains for us is to show that the terms contributing in the sums on the right side 
of ( ]C.12I) and corresponding to the case Dn-, d > 2m — k and d + s < k — 2r, fulfill the 



annihilation thesis, too. 

It's reasonable to call the sum over all terms on the right side of (1C.12I) which fulfill 
the condition for the case Dn as D^2ms^ following calculations, because the total 

propagation time is for the term -Dfj*2m s'* ^^^^^ t = {n + l)At, but now this term is 
related to a case differentiation done for a calculation where we wrote out the operator 
W{n) explicitly. We can note these terms -D^^^gm s'' "with the conditions d > 2m — k and 
d + s < k — 2r and (IC.12I) in the following way, where we regard that the inequations 
s<k — 2r — d<k — 2r — {2m — k) = 2{k — r — m) hold, and moreover, we substitute 
for systematical reasons the sum index r by r„_|_i, because 2r is the power in which the 
operator W{n + 1) appears in (]C.12I) : 



k-m min[2{k-rn+i-m), 2(m-r„+i)] 

^Sr^v? = i-ir^t^- E w{n + ir-^^ mod 2, >< 

rn+l=0 s=0 
min[2{rn—rn+i)—s, fc— 2r„^i, fc— 2r„+i— s] 



X 



X 



d=max{2m—k, 2m—k—s, 0) 
n-1 



q=l 
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n-2 



7'^{n-l,2(m-r„+i)-d-s) P-0 



(n — l,2(m — rj^_j_]^) — d — s) 



iHoAtj 



#(0) 



(C.17) 



Now we consider like in the calculation of ( 1C.9|) that for the oscillatory orders holds 
k < 2m < 2k, thus, we can simplify the max and min functions in the sum limits over s 
and d in (lC.17p . and by taking (|3]) into account and renaming 2s by r„ and dhj j, we 
get as a result for 



D 



(tn + l,fc) 

n,2m,S 



2m,S ■ 

k-m k—rn+i—m 
rn+l=0 rn=0 
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Now we compare flC.9P to (IC.ISP and realize that we can iterate the implication from 
(IC.9P to f lC.lSP and thus, get with the abbreviation 
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a general expression for -D^*^!^^^ for a = 2, 3 
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The idea to prove the annihilation thesis is now that if we can prove that one term 
D^a2m s^ for any value for a = 2, 3, . . . , + 1 is proportional to a global sign factor 
(— we have accomplished the proof of the annihilation thesis. The reason for this 
is that because of the results attained in the case differentiations before, which we did in 
the implication from the oscillatory order A^^^^^ 2m ^ < 2m < 2k to the term D^2m S^ ^ 
the term A^^+i_2m,s is for a proportionality of the term D^^^mS^ the global sign factor 
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, too. 



For this task we calculate with ( !Cl20l) the term D^^"^^'^^- 
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fc-2 S2 

j=2m— A; 
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Moreover, since it is evident that the equations 
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hold, thus, we can derive for D^^^mp ^^^^ 
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Then with (]B.7P we can carry out the sum over j and get as a final result for -D2 2m 5 • 
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\A:— m 



We perceive from f lC.24p that 2m's proportional to a global sign factor ( 
thus, as we discussed before, that means that we have proven for the oscillatory orders 
N^2mSi^ < 2m < 2k the annihilation thesis that these orders are proportional to a 
global sign factor (— and that all terms with an sign factor (—i^'^-^+i disappear. 
□ 
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